/*==========================================================================================     * Date: Sept 2018
	* Last Modified: 23/05/2017
	* Purpose: 
		- Border specification with weights/nweights and discontinuity graph
		- creates function GetDistDummy(string var, string all, string social) to create  function generating total weighted employment within radius per point in the map GetDistDummy(labour variable, Total obs var , social group var)
		- creates function GetDistTotal(string all, string social) to create function generating total weighted observations within radius per point in the map GetDistDummy(Total obs var , social group var)

	* DATASETS Used:
		- "$NAPP/MolaModified_NAPP_JAG4new"
		
	* DOFILES used:
		- "$CODE/bordersocial.do"
				
	* DATASETS Created:
		- "$NAPP/Proof2new.dta"
		- "$NAPP/Proof2Rnew.dta"
==========================================================================================*/


gl CODE "/Users/`c(username)'/Dropbox/NetworkParish2016/Occupational_Choice/Submission/Final_ReStat/Codes/C_ANALYSIS"
*gl CODE "$ROOT/NetworkParish2016/Occupational_Choice/Codes/Stata Code/FINAL/HPC_CODE"
gl NAPP "/Users/`c(username)'/Dropbox/NetworkParish2016/Occupational_Choice/Submission/Final_ReStat/Codes/NAPP"
*gl NAPP "$ROOT/NetworkParish2016/Occupational_Choice/Data/NAPP_address"
gl RESULTS "/Users/`c(username)'/Dropbox/NetworkParish2016/Occupational_Choice/Submission/Final_ReStat/Codes/Results"
*gl RESULTS "$ROOT/NetworkParish2016/Occupational_Choice/Results/EEA"
gl MOLA "/Users/`c(username)'/Dropbox/NetworkParish2016/Occupational_Choice/Submission/Final_ReStat/Codes/MOLA"

clear all
cap log close
set more off
set matsize 11000
set maxvar 8000
set emptycells drop

/*=============================*/
/* Create sample for the study */
/*=============================*/

use "$NAPP/MolaModified_NAPP_JAG4new", clear

do "$CODE/bordersocial.do"

levelsof ID_socialborder if bordersocial==1, local(dropborder)

drop if bordersocial==1 //I eliminate this because they do not have correct measures

*Some definitions
sort HHNUMhead PERNUMhead
gen child = (tot_child>=1)
gen nchild=tot_child 
recode SERVANGB (99=.)
gen nservant = SERVANGB

recode labforce (1=1) (2=0) (9=.)
recode SEX (1=1) (2=0) (9=.)

* Sample restrictions
* Keep sample of working age individual not living in an institution
keep if AGE>=15 & AGE<=60
rename employednew labforcenew
drop if labforcenew==.
drop if SEX==.
drop if MARST==9
drop if nativity==9
*To keep only household head, if we want to have everybody we should eliminate this
keep if hhhead==1
*drop if they live in an institution
gen residential=(INSTITGB==1 | INSTITGB==5)
drop if residential==0
*drop if they belong to the servants hh (they also live where they work)
drop if HHservantunknown==1

rename objectid_1 OBJECTID_1
cap drop _merge

gen married = .
replace married = 1 if MARST==1 | MARST==2 | MARST==4
replace married = 0 if MARST==6

gen foreign = .
replace foreign = 1 if nativity==3 | nativity==2 | nativity==8
replace foreign = 0 if nativity==1
*drop foreigners
drop if foreign==1
gen native = (foreign==0)

gen stayerp = (migrant==1)
gen leaverp = (stayerp==0)

gen stayerc = (migrant<=2)
gen leaverc = (stayerc==0)

recode UNMARDAU (99=.) 
recode UNMARSON (99=.)
recode MARRYDAU (99=.)
recode MARRYSON (99=.)
recode UNMARKID (99=.)

gen single_f = (married==0 & SEX==2)
gen single_m = (married==0 & SEX==1)

gen age15_16 = (AGE>=15 & AGE<=16)
gen age17_20 = (AGE>=17 & AGE<=20)
gen age21_60 = (AGE>=21 & AGE<=60)

gen women = (SEX==0)
gen men = (SEX==1)

gen one=1
local subgroups "one age15_16 age17_20 age21_60"
local labours "labforce class1 class2 class3 class5 labforcenew classnew0 classnew1 classnew2 classnew3 classnew4 classnew5 classnew6 classnew7 AGE SEX married nchild nservant stayerp stayerc"
foreach subgroup of local subgroups {
	foreach labour of local labours {
		gen `labour'_`subgroup'= `labour'*`subgroup'
	}
}

local subgroups "women men"
local labours "labforce class1 class2 class3 class5 labforcenew classnew0 classnew1 classnew2 classnew3 classnew4 classnew5 classnew6 classnew7 AGE married nchild nservant stayerp stayerc"
foreach subgroup of local subgroups {
	foreach labour of local labours {
		gen `labour'_`subgroup'= `labour'*`subgroup'
	}
}

sort OBJECTID_1
cap gen one=1
collapse (sum) one labforce* class1* class2* class3* class5* classnew0* classnew1* classnew2* classnew3* classnew4* classnew5* classnew6* classnew7* ///
women men AGE* SEX* married* nchild* nservant* foreign* stayer* age15_16* age17_20* age21_60*  ///
(mean) social civil Buffer100 Buffer80 Buffer60 Buffer50 Buffer40 Buffer30 ID_socialborder IDSocial_L IDSocial_R, by(OBJECTID_1)
sort OBJECTID_1
tempfile Reduced
save `Reduced', replace
save "$MOLA\ReducedSamplenew", replace

**** TO DO: keep objectid_1 longitude latitude
**** TemporaryPointsCoordinatesv1.do

/*===========================================*/
/* Create distance matrix to weight distance */
/*===========================================*/

use "$MOLA\points1881DistMatrix.dta", clear
sort OBJECTID_1
merge 1:1 OBJECTID_1 using "$MOLA\ReducedSamplenew"
levelsof OBJECTID_1 if _merge~=3, local(nomatched)
egen NOOK=anymatch(OBJECTID_1), values(`nomatched')
drop if NOOK
cou
foreach nomatch of local nomatched {
	cap drop d`nomatch' 
}
cap drop _merge NOOK
sum OBJECTID_1
local idmin=r(min)
local idmax=r(max)
order OBJECTID_1 d`idmin'-d`idmax'


**To create  Function generating total weighted employment within radius per point in the map GetDistDummy(labour variable, Total obs var , social group var)
set more off
cap mata mata drop GetDistDummy()
cap mata mata clear
mata :
function GetDistDummy(string var, string all, string social) 
{
nobs = st_nobs()
nobs1=nobs+1
st_view(Dist=.,.,2..nobs1)
st_view(Out=.,.,var)
st_view(All=.,.,all)
st_view(Social=.,.,social)
SSocial=J(nobs,1,1)#Social'
DSocial=Social:==SSocial
DDist=.
Disttemp=Dist:^(-1)
Disttemp=editmissing(Disttemp,0)
MaxDisttemp=rowmax(Disttemp)
Diag=diag(MaxDisttemp)
DDist=Disttemp+Diag
DDistSocial=DDist:*DSocial
DDistlab=.
DDistlab= DDist * Out
DDistsoclab= DDistSocial * Out
DDistdiflab=DDistlab-DDistsoclab
Distvar=st_addvar("double","T"+var)
st_store((1,rows(DDistlab)),Distvar,DDistlab)
Distsocvar=st_addvar("double","T"+social+var)
st_store((1,rows(DDistsoclab)),Distsocvar,DDistsoclab)
Distdifvar=st_addvar("double","TN"+social+var)
st_store((1,rows(DDistdiflab)),Distdifvar,DDistdiflab)
DDistall= DDist * All
DDistsocall= DDistSocial * All
DDistdifall=DDistall-DDistsocall
DDistmean=DDistlab:/DDistall
Distmeanvar=st_addvar("double","M"+var)
st_store((1,rows(DDistmean)),Distmeanvar,DDistmean)
DDistsocmean=DDistsoclab:/DDistsocall
DDistdifmean=DDistdiflab:/DDistdifall
Distsocmeanvar=st_addvar("double","M"+social+var)
st_store((1,rows(DDistsocmean)),Distsocmeanvar,DDistsocmean)
Distdifmeanvar=st_addvar("double","MN"+social+var)
st_store((1,rows(DDistdifmean)),Distdifmeanvar,DDistdifmean)
}
end

**To create  Function generating total weighted observations within radius per point in the map GetDistDummy(Total obs var , social group var)
cap mata mata drop GetDistTotal()
mata :
function GetDistTotal(string all, string social) 
{
nobs = st_nobs()
nobs1=nobs+1
st_view(Dist=.,.,2..nobs1)
st_view(All=.,.,all)
st_view(Social=.,.,social)
SSocial=J(nobs,1,1)#Social'
DSocial=Social:==SSocial
DDist=.
Disttemp=Dist:^(-1)
Disttemp=editmissing(Disttemp,0)
MaxDisttemp=rowmax(Disttemp)
Diag=diag(MaxDisttemp)
DDist=Disttemp+Diag
DDistSocial=DDist:*DSocial
DDistall= DDist * All
Distallvar=st_addvar("double","T"+all)
st_store((1,rows(DDistall)),Distallvar,DDistall)
DDistsocall= DDistSocial * All
DDistdifall=DDistall-DDistsocall
Distsocallvar=st_addvar("double","T"+social+all)
st_store((1,rows(DDistsocall)),Distsocallvar,DDistsocall)
Distdifallvar=st_addvar("double","TN"+social+all)
st_store((1,rows(DDistdifall)),Distdifallvar,DDistdifall)
DistWeightvar=st_addvar("double","W"+social+all)
st_store((1,rows(MaxDisttemp)),DistWeightvar,MaxDisttemp)
}
end

set more off
**maximum distance within same social 1724.266113
local subgroups "one age15_16 age17_20 age21_60"
foreach subgroup of local subgroups {
	foreach var of varlist labforce class1 class2 class3 class5 labforcenew classnew0 classnew1 classnew2 classnew3 classnew4 classnew5 classnew6 classnew7 AGE SEX married nchild nservant stayerp stayerc {
		mata GetDistDummy("`var'_`subgroup'","`subgroup'","social")
	}
	mata GetDistTotal("`subgroup'","social")	
}

local subgroups "women men"
foreach subgroup of local subgroups {
	foreach var of varlist labforce class1 class2 class3 class5 labforcenew classnew0 classnew1 classnew2 classnew3 classnew4 classnew5 classnew6 classnew7 AGE married nchild nservant stayerp stayerc {
		mata GetDistDummy("`var'_`subgroup'","`subgroup'","social")
	}
	mata GetDistTotal("`subgroup'","social")	
}

drop d`idmin'-d`idmax'	
keep OBJECTID_1 T* M* Wsocial* social civil Buffer100 Buffer80 Buffer60 Buffer50 Buffer40 Buffer30 ID_socialborder IDSocial_L IDSocial_R
sort OBJECTID_1
gen objectid_1=OBJECTID_1
save "$NAPP/WeightedVariablesnew.dta", replace

/*=======================================================*/
/* Create variables with our dataset */
/*=======================================================*/

set more off
use "$NAPP/WeightedVariablesnew.dta", clear
merge 1:n objectid_1 using "$NAPP/MolaModified_NAPP_JAG4new", keep(matched) generate(newmerge)
rename employednew labforcenew
save "$NAPP/Proofnew.dta", replace


set more off
use "$NAPP/Proofnew.dta", clear

// Sample 

recode labforce (1=1) (2=0) (9=.)
recode SEX (1=1) (2=0) (9=.)

sort HHNUMhead PERNUMhead
gen child = (tot_child>=1)
gen nchild=tot_child 

recode SERVANGB (99=.)
gen nservant = SERVANGB

* Keep sample of working age individual not living in an institution
keep if AGE>=15 & AGE<=60
drop if labforcenew==.
drop if SEX==.
drop if MARST==9
drop if nativity==9
*To keep only household head, if we want to have everybody we should eliminate this
keep if hhhead==1
*drop if they live in an institution
gen residential=(INSTITGB==1 | INSTITGB==5)
drop if residential==0
*drop if they belong to the servants hh (they also live where they work)
drop if HHservantunknown==1

* Create some explanatory variables
gen married = .
replace married = 1 if MARST==1 | MARST==2 | MARST==4
replace married = 0 if MARST==6

gen foreign = .
replace foreign = 1 if nativity==3 | nativity==2 | nativity==8
replace foreign = 0 if nativity==1
drop if foreign==1
gen native = (foreign==0)

gen stayerp = (migrant==1)
gen leaverp = (stayerp==0)

gen stayerc = (migrant<=2)
gen leaverc = (stayerc==0)

recode UNMARDAU (99=.) 
recode UNMARSON (99=.)
recode MARRYDAU (99=.)
recode MARRYSON (99=.)
recode UNMARKID (99=.)

gen single_f = (married==0 & SEX==2)
gen single_m = (married==0 & SEX==1)

cap gen AGE2 = AGE^2

// In same occupation as father

// In same occupation as a family member

// AGE AGE2 married nchild SERVANGB foreign stayer
cap gen one=1
cap gen id = 1

// A1 - Local employment rate
bys social civil: egen tot_empl_s = sum(labforce)
bys social civil: egen tot_emplnew_s = sum(labforcenew)
bys social civil: egen tot_people_s = sum(id)

// A2 - Employment for all except members of the household
bys HHNUMhead: egen empl_hh = sum(labforce)
bys HHNUMhead: egen emplnew_hh = sum(labforcenew)
gen tot_empl_hh = tot_empl_s - empl_hh
gen tot_emplnew_hh = tot_emplnew_s - emplnew_hh

bys HHNUMhead: egen people_hh = sum(id)
gen tot_people_hh = tot_people_s - people_hh

* Non weighted share in vars w/o hh obs
foreach var in AGE SEX married nchild nservant foreign native stayerp stayerc leaverp leaverc {
	bys civil social: egen tot_`var' = sum(`var')
	bys HHNUMhead: egen tot_`var'_hh = sum(`var')
	gen mean_nw_`var'= (tot_`var'-tot_`var'_hh)/tot_people_hh
}

* Non weighted share in labforce w/o hh obs
cap gen sh_nw_labforce_hh = tot_empl_hh / tot_people_hh
cap label variable sh_nw_labforce_hh "Share nw empl"

cap gen sh_nw_labforcenew_hh = tot_emplnew_hh / tot_people_hh
cap label variable sh_nw_labforcenew_hh "Share nw empl new"

// A3 - Employment for men only except all members of the household
gen men = (SEX==1)
bys social civil: egen tot_empl_men = sum(labforce*men)
bys social civil: egen tot_emplnew_men = sum(labforcenew*men)
bys social civil: egen tot_men = sum(men)

bys HHNUMhead: egen empl_men_hh = sum(labforce*men)
bys HHNUMhead: egen emplnew_men_hh = sum(labforcenew*men)
gen tot_empl_men_hh = tot_empl_men - empl_men_hh
gen tot_emplnew_men_hh = tot_emplnew_men - emplnew_men_hh

bys HHNUMhead: egen men_hh = sum(men)
gen tot_men_hh = tot_men - men_hh

cap gen sh_nw_labforce_men_hh = tot_empl_men_hh / tot_men_hh
label variable sh_nw_labforce_men_hh "Share nw empl men"
cap gen sh_nw_labforcenew_men_hh = tot_emplnew_men_hh / tot_men_hh
label variable sh_nw_labforcenew_men_hh "Share nw empl new men"

// A4 - Employment for women only except all members of the household
gen women = (SEX==0)
bys social civil: egen tot_empl_women = sum(labforce*women)
bys social civil: egen tot_emplnew_women = sum(labforcenew*women)
bys social civil: egen tot_women = sum(women)

bys HHNUMhead: egen empl_women_hh = sum(labforce*women)
bys HHNUMhead: egen emplnew_women_hh = sum(labforcenew*women)
gen tot_empl_women_hh = tot_empl_women - empl_women_hh
gen tot_emplnew_women_hh = tot_emplnew_women - emplnew_women_hh

bys HHNUMhead: egen women_hh = sum(women)
gen tot_women_hh = tot_women - women_hh

gen sh_nw_labforce_women_hh = tot_empl_women_hh / tot_women_hh
label variable sh_nw_labforce_women_hh "Share nw empl women"
gen sh_nw_labforcenew_women_hh = tot_emplnew_women_hh / tot_women_hh
label variable sh_nw_labforcenew_women_hh "Share nw empl new women"

// B1 - Share of employed in your own parish
levelsof social, local(social_parish)
foreach y of local social_parish {
      gen x_`y' = (social==`y')
      egen tot_empl_`y' = sum(x_`y'*labforce) 
	  egen tot_emplnew_`y' = sum(x_`y'*labforcenew) 
      egen tot_pop_`y' = sum(x_`y')
      gen sh_empl_`y' = tot_empl_`y' / tot_pop_`y'
	  gen sh_emplnew_`y' = tot_emplnew_`y' / tot_pop_`y'
      drop tot_empl_`y' tot_emplnew_`y' tot_pop_`y' x_`y' 
}

// B2 - Share of employed in control parishes (not your own parish)
gen control_parish = IDSocial_L if IDSocial_L!=social
replace control_parish = IDSocial_R if IDSocial_R!=social

gen sh_nw_labforce_control = .
gen sh_nw_labforcenew_control = .

foreach y of local social_parish{
      replace sh_nw_labforce_control = sh_empl_`y' if control_parish==`y'
	  replace sh_nw_labforcenew_control = sh_emplnew_`y' if control_parish==`y'
}

// B3 - Share of people in each class by social parish non-weighted
foreach occ in class1 class2 class3 class5 {
	bys social civil: egen tot_`occ' = sum(`occ')
	bys HHNUMhead: egen `occ'_hh = sum(`occ'*labforce)
	replace tot_`occ' = tot_`occ' - `occ'_hh
	gen sh_nw_`occ'_hh = tot_`occ' /  tot_people_hh
}

foreach occ in classnew0 classnew1 classnew2 classnew3 classnew4 classnew5 classnew6 classnew7 {
	bys social civil: egen tot_`occ' = sum(`occ')
	bys HHNUMhead: egen `occ'_hh = sum(`occ'*labforcenew)
	replace tot_`occ' = tot_`occ' - `occ'_hh
	gen sh_nw_`occ'_hh = tot_`occ' /  tot_people_hh
}
// Share of people in each age group

gen age15_16 = (AGE>=15 & AGE<=16)
gen age17_20 = (AGE>=17 & AGE<=20)
gen age21_60 = (AGE>=21 & AGE<=60)

foreach x in 15_16 17_20 21_60 {
	bys social civil: egen tot_age`x' = sum(age`x')
	bys social civil: egen tot_empl_age`x' = sum(labforce*age`x')
	bys social civil: egen tot_emplnew_age`x' = sum(labforcenew*age`x')
    
	bys HHNUMhead: egen age`x'_hh = sum(age`x')
	bys HHNUMhead: egen empl_age`x'_hh = sum(age`x'*labforce)
	bys HHNUMhead: egen emplnew_age`x'_hh = sum(age`x'*labforcenew)
      
	replace tot_age`x' = tot_age`x' - age`x'_hh
	replace tot_empl_age`x' = tot_empl_age`x' - empl_age`x'_hh
	replace tot_emplnew_age`x' = tot_emplnew_age`x' - emplnew_age`x'_hh

	gen sh_nw_labforce_age`x'_hh = tot_empl_age`x' /  tot_age`x'
	gen sh_nw_labforcenew_age`x'_hh = tot_emplnew_age`x' /  tot_age`x'
}

// Share of people in each occupation inside the family

// Share of men in each occupation inside the family

// C - Distance weighted shares 

set more off
local subgroups "one women men age15_16 age17_20 age21_60"

foreach subgroup of local subgroups {
	cap drop temphh`subgroup'
	cap drop temppoint`subgroup'
	bys HHNUMhead: egen double temphh`subgroup'= sum(`subgroup'*Wsocial`subgroup')
	bys objectid_1: egen double temppoint`subgroup'= sum(`subgroup'*Wsocial`subgroup')
	
	foreach var of varlist labforce class1 class2 class3 class5 labforcenew classnew0 classnew1 classnew2 classnew3 classnew4 classnew5 classnew6 classnew7 AGE married nchild nservant stayerp stayerc {
		cap drop temphh`var'`subgroup'
		cap drop temppoint`var'`subgroup'
		bys HHNUMhead: egen double temphh`var'`subgroup'= sum(`var'*`subgroup'*Wsocial`subgroup')
		bys objectid_1: egen double temppoint`var'`subgroup'= sum(`var'*`subgroup'*Wsocial`subgroup')

		gen double Totsoc`var'`subgroup'hh=Tsocial`var'_`subgroup'-temphh`var'`subgroup'
		replace Totsoc`var'`subgroup'hh=round(Totsoc`var'`subgroup'hh*100000)/100000 
		cap gen double Totsoc`subgroup'hh=Tsocial`subgroup'-temphh`subgroup'
		replace Totsoc`subgroup'hh=round(Totsoc`subgroup'hh*100000)/100000 
		gen double Measoc`var'`subgroup'hh=Totsoc`var'`subgroup'hh/Totsoc`subgroup'hh
		gen double MeaNsoc`var'`subgroup'hh=MNsocial`var'_`subgroup'
		
		gen double Totsoc`var'`subgroup'pp=Tsocial`var'_`subgroup'-temppoint`var'`subgroup'
		replace Totsoc`var'`subgroup'pp=round(Totsoc`var'`subgroup'pp*100000)/100000 
		cap gen double Totsoc`subgroup'pp=Tsocial`subgroup'-temppoint`subgroup'
		replace Totsoc`subgroup'pp=round(Totsoc`subgroup'pp*100000)/100000 
		gen double Measoc`var'`subgroup'pp=Totsoc`var'`subgroup'pp/Totsoc`subgroup'pp
		gen double MeaNsoc`var'`subgroup'pp=MNsocial`var'_`subgroup'
	}
}

local subgroups "one age15_16 age17_20 age21_60"

foreach subgroup of local subgroups {
	cap drop temphh`subgroup'
	cap drop temppoint`subgroup'
	bys HHNUMhead: egen double temphh`subgroup'= sum(`subgroup'*Wsocial`subgroup')
	bys objectid_1: egen double temppoint`subgroup'= sum(`subgroup'*Wsocial`subgroup')
	foreach var of varlist SEX {
		cap drop temphh`var'`subgroup'
		cap drop temppoint`var'`subgroup'
		bys HHNUMhead: egen double temphh`var'`subgroup'= sum(`var'*`subgroup'*Wsocial`subgroup')
		bys objectid_1: egen double temppoint`var'`subgroup'= sum(`var'*`subgroup'*Wsocial`subgroup')

		gen double Totsoc`var'`subgroup'hh=Tsocial`var'_`subgroup'-temphh`var'`subgroup'
		replace Totsoc`var'`subgroup'hh=round(Totsoc`var'`subgroup'hh*100000)/100000 
		cap gen double Totsoc`subgroup'hh=Tsocial`subgroup'-temphh`subgroup'
		replace Totsoc`subgroup'hh=round(Totsoc`subgroup'hh*100000)/100000 
		gen double Measoc`var'`subgroup'hh=Totsoc`var'`subgroup'hh/Totsoc`subgroup'hh
		gen double MeaNsoc`var'`subgroup'hh=MNsocial`var'_`subgroup'
		
		gen double Totsoc`var'`subgroup'pp=Tsocial`var'_`subgroup'-temppoint`var'`subgroup'
		replace Totsoc`var'`subgroup'pp=round(Totsoc`var'`subgroup'pp*100000)/100000 
		cap gen double Totsoc`subgroup'pp=Tsocial`subgroup'-temppoint`subgroup'
		replace Totsoc`subgroup'pp=round(Totsoc`subgroup'pp*100000)/100000 
		gen double Measoc`var'`subgroup'pp=Totsoc`var'`subgroup'pp/Totsoc`subgroup'pp
		gen double MeaNsoc`var'`subgroup'pp=MNsocial`var'_`subgroup'
	}
}

// 

cap drop temp* sh_empl_* sh_emplnew_*
save "$NAPP/Proof1new.dta", replace
*
set more off
use "$NAPP/Proof1new.dta", clear

set more off
foreach subgroup in women men age15_16 age17_20 age21_60 {
	foreach var of varlist labforce class1 class2 class3 class5 labforcenew classnew0 classnew1 classnew2 classnew3 classnew4 classnew5 classnew6 classnew7 {
	rename Measoc`var'`subgroup'hh sh_`var'_`subgroup'_hh
	}
	foreach var of varlist AGE married nchild nservant stayerp stayerc {
	rename Measoc`var'`subgroup'hh mean_`var'_`subgroup'
	}
}

foreach subgroup in age15_16 age17_20 age21_60 {
	foreach var of varlist SEX {
	rename Measoc`var'`subgroup'hh mean_`var'_`subgroup'
	}
}

foreach subgroup in one {
	foreach var of varlist labforce class1 class2 class3 class5 labforcenew classnew0 classnew1 classnew2 classnew3 classnew4 classnew5 classnew6 classnew7 {
	cap drop sh_`var'_hh
	rename Measoc`var'`subgroup'hh sh_`var'_hh
	}
	foreach var of varlist AGE SEX married nchild nservant stayerp stayerc {
	cap rename mean_`var' mean_nw_`var'
	rename Measoc`var'`subgroup'hh mean_`var'
	}
}

*To create control parish
levelsof social, local(social_parish)
foreach y of local social_parish {
	gen x_`y' = 1 if social==`y'
	foreach var of varlist labforce class1 class2 class3 class5 labforcenew classnew0 classnew1 classnew2 classnew3 classnew4 classnew5 classnew6 classnew7 {
		egen sh_`var'_hh_`y' = mean(x_`y'*sh_`var'_hh)
	}
	drop x_`y'
}

set more off
levelsof social, local(social_parish)
foreach var of varlist labforce class1 class2 class3 class5 labforcenew classnew0 classnew1 classnew2 classnew3 classnew4 classnew5 classnew6 classnew7 {
	gen sh_`var'_control_hh=.
}
foreach y of local social_parish {
	foreach var of varlist labforce class1 class2 class3 class5 labforcenew classnew0 classnew1 classnew2 classnew3 classnew4 classnew5 classnew6 classnew7 {
		replace sh_`var'_control_hh = sh_`var'_hh_`y' if control_parish==`y'
		drop sh_`var'_hh_`y'
	}
}
foreach var of varlist labforce class1 class2 class3 class5 labforcenew classnew0 classnew1 classnew2 classnew3 classnew4 classnew5 classnew6 classnew7 {
	gen dgood`var'=(sh_`var'_hh>sh_`var'_control_hh)
	replace dgood`var'=. if sh_`var'_control_hh==.
}
// Labels
set more off
label variable AGE "age"
label variable AGE2 "age2"
label variable SEX "male"
label variable married "married"
label variable nchild "n children"
label variable nservant "n servants"
label variable SERVANGB "servants"
label variable foreign "foreigner"
label variable stayerp "stayer in parish"
label variable stayerc "stayer in county"
*label variable sh_empl_control "share of empl in control"
cap label variable sh_class1 "sh professional"
cap label variable sh_class2 "sh domestic"
cap label variable sh_class3 "sh commercial"
*cap label variable sh_class4 "sh agricultural"
cap label variable sh_class5 "sh industrial"

cap label variable sh_classnew0 "sh unemployed"
cap label variable sh_classnew1 "sh professional"
cap label variable sh_classnew2 "sh domestic"
cap label variable sh_classnew3 "sh commercial"
cap label variable sh_classnew4 "sh ind-artisan"
cap label variable sh_classnew5 "sh ind-builder"
cap label variable sh_classnew6 "sh ind-food"
cap label variable sh_classnew7 "sh ind-services"

cap label variable sh_age15_16 "sh age 15-16"
cap label variable sh_age17_20 "sh age 17-20"
cap label variable sh_age21_60 "sh age 21-60"


// Drop borders ID from extreme borders (XXX We should comment dropborder if we run the full code XXX)
cap egen dropboundaries = anymatch(ID_socialborder), values(`dropborder')
local dropborder="1 2 3 4 5 6 7 8 9 10 12 13 15 17 18 20 23 25 26 28 30 31 33 34 35 36 38 39 40 41 42 43 44 45 46 47 48 53 54 58 59 60 61 62 63 67 71 72 77 78 79 82 84 91 92 97 101 104 105 108 109 110 121 138 141 148 158 171 172 184 214 216 225 247 252 253 272 305 314 319 333 336 344 369 382 463 464 521 531 568 573 582 583 665 682 683 685 754 755 775 796 804 817 827 831 832 859 863 864 865 867 883 887 895 901 902 934 950 955 957 960 962 964 970 991 992 993 994 995 998 1003 1004 1010 1011 1012 1013 1014 1020 1023 1024 1025 1026 1027 1029 1031 1033 1036 1041 1042 1043 1044 1045 1047 1049 1051 1052 1054 1055 1056 1058 1059 1060 1061 1062 1063 1064 1066 1068 1069 1070 1071 1073 1074 1075 1077 1078 1079 1080 1081 1082 1084 1085 1086 1087 1088 1089 1090 1091 1093 1094 1095 1096 1097 1098"
cap egen dropboundaries = anymatch(ID_socialborder) , values(`dropborder')
drop if dropboundaries==1

// Taking away civil parish== social parish borders

bys ID_socialborder civil: gen problem = _n==1
bys ID_socialborder: egen problem2 = sum(problem)
*drop if problem2>1

// Taking away borders with only one social parish (therefore no control)

bys ID_socialborder social: gen problemsocial = _n==1
bys ID_socialborder: egen problemsocial2 = sum(problemsocial)
*drop if problemsocial2==1

*Discontinuity graph
set more off
gen minusdist_socialborder=-1*dist_socialborder

**** TO DO: keep objectid_1 longitude latitude
**** TemporaryPointsCoordinatesv1.do

save "$NAPP/Proof2new.dta", replace

keep HHNUMhead PERNUMhead objectid_1 sh_*_hh mean_* nchild stayerp stayerc nservant SEX married AGE problem2 problemsocial2 Buffer100 Buffer80 Buffer60 Buffer50 Buffer40 classgb classnew civil social ID_socialborder
*save in version 12 so R canimport it
saveold "$NAPP/Proof2Rnew.dta", replace version(12)